library(fixest) # to demean variables
library(multiwayvcov) # to cluster SEs
library(sandwich) # to get HC2 SEs
library(stargazer) # to generate tables
library(mvtnorm) # to simulate betas and vcov matrix

rm(list = ls())

# load data
load(file = 'chData.RData')

# Put each dataset in a different object
ch_reform <- ch_list[[1]]
ch_refTerm <- ch_list[[2]]
ch_pre_post <- ch_list[[3]]
ch_refTermMChanges <- ch_list[[4]]
ch_refTermParty <- ch_list[[5]]
ch_refTermMChangesParty <- ch_list[[6]]

################################# 
####### Models using M ##########
#################################


# C stands for congress
# G stands for gender
# M stands for magnitude

# Model 1: G*M C1,C2,C3 (all data)
data_m1 <- demean(women_pct + magAnalysis + female + magAnalysisWoman
	~ session, data = ch_refTerm)
data_m1 <- cbind(data_m1, bill_author = ch_refTerm$bill_author)
m1 <- lm(women_pct ~ magAnalysis + female + magAnalysisWoman - 1, data = data_m1)
vcov1 <- cluster.vcov(m1, cluster = data_m1$bill_author)
s1  <- sqrt(diag(vcov1))

# Model 2: G*M C1,C2,C3 (M changes at the middle of C2)
data_m2 <- demean(women_pct + magAnalysis + female + magAnalysisWoman
	~ session + reform, data = ch_refTermMChanges)
data_m2 <- cbind(data_m2, bill_author = ch_refTermMChanges$bill_author)
m2 <- lm(women_pct ~ magAnalysis + female + magAnalysisWoman - 1, data = data_m2)
vcov2 <- cluster.vcov(m2, cluster = data_m2$bill_author)
s2  <- sqrt(diag(vcov2))


# Model3: G*M C1,C3
data_m3 <- demean(women_pct + magAnalysis + female + magAnalysisWoman
	~ session, data = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022')))
data_m3 <- cbind(data_m3, bill_author = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022'))$bill_author)
m3 <- lm(women_pct ~ magAnalysis + female + magAnalysisWoman - 1, data = data_m3)
vcov3 <- cluster.vcov(m3, cluster = data_m3$bill_author)
s3  <- sqrt(diag(vcov3))

# Model 4: C2
inSample <- names(table(ch_pre_post$bill_author)[table(ch_pre_post$bill_author) == 2])
ch_pre_post <- subset(ch_pre_post, bill_author %in% inSample)
data_m4 <- fixest::demean(women_pct + magAnalysis + female + magAnalysisWoman
	~ bill_author, data = ch_pre_post)
data_m4 <- cbind(data_m4, bill_author = ch_pre_post$bill_author)
m4 <- lm(women_pct ~ magAnalysis + magAnalysisWoman - 1, data = data_m4)
vcov4 <- cluster.vcov(m4, cluster = data_m4$bill_author)
s4  <- sqrt(diag(vcov4))

# Model 5: G*M C1,C2,C3 (only repeated legislators)
data_m5 <- demean(women_pct + magAnalysis + magAnalysisWoman
	~ bill_author + session, data = subset(ch_refTerm, allSessions == 1))
data_m5 <- cbind(data_m5, bill_author = subset(ch_refTerm, allSessions == 1)$bill_author)
m5 <- lm(women_pct ~ magAnalysis + magAnalysisWoman - 1, data = data_m5)
vcov5 <- cluster.vcov(m5, cluster = data_m5$bill_author)
s5  <- sqrt(diag(vcov5))

# Model 6: G*M C1, C2, C3 (M changes at the middle of C2) (only repeated legislators)
data_m6 <- demean(women_pct + magAnalysis + magAnalysisWoman
	~ session + reform + bill_author, data = subset(ch_refTermMChanges, allSessions == 1))
data_m6 <- cbind(data_m6, bill_author = subset(ch_refTermMChanges, allSessions == 1)$bill_author)
m6 <- lm(women_pct ~ magAnalysis + magAnalysisWoman - 1, data = data_m6)
vcov6 <- cluster.vcov(m6, cluster = data_m6$bill_author)
s6  <- sqrt(diag(vcov6))

# Model 7: G*M C1,C3 (only repeated legislators)
data_m7 <- demean(women_pct + magAnalysis + magAnalysisWoman
	~ bill_author + session, data = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022') & l10_l18 == 1))
data_m7 <- cbind(data_m7, bill_author = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022') & l10_l18 == 1)$bill_author)
m7 <- lm(women_pct ~ magAnalysis + magAnalysisWoman - 1, data = data_m7)
vcov7 <- cluster.vcov(m7, cluster = data_m7$bill_author)
s7  <- sqrt(diag(vcov7))

# Model 8: C3
data_m8 <- subset(ch_refTerm, session == '2018-2022')
m8 <- lm(women_pct ~ magAnalysis + female + magAnalysisWoman, data = data_m8)
vcov8 <- vcovHC(m8, 'HC2')
s8  <- sqrt(diag(vcov8))

# Generate Table C.9
stargazer(m1, m2, m3, m4, m5, m6, m7, m8, 
	se = list(s1, s2, s3, s4, s5, s6, s7, s8), 
	covariate.labels = c('M', 'Woman', 'M x Woman'), 
	add.lines = list(c('FE by Legislative Term', 'Yes', 'Yes', 'Yes', 'No', 'Yes', 'Yes', 'Yes', 'No'),
		c('FE by Reform', 'No', 'Yes', 'No', 'No', 'No', 'Yes', 'No', 'No'),
		c('FE by Legislator', 'No', 'No', 'No', 'Yes', 'Yes', 'Yes', 'Yes', 'No')),
	dep.var.labels = "Bills' Portfolio on Women's Issues (%)",
	omit.stat = c("ser",'f'),
	star.cutoffs = c(0.10, 0.05, 0.01),
	no.space = TRUE,
	single.row = FALSE,	
	notes.append = FALSE,
	title = "Association Between Legislative Portfolio, Gender, and District Magnitude–2014-2021 Cámara de Diputadas y Diputados de Chile",
	label = 't:mag'
)

#################### Substantive Effect ########################
## Scenarios: Women, M = 2
W2 <- data.frame(mag = 2, female = 1, magAnalysisWoman = 2 * 1)
## Scenarios: Women, M = 3 (post-reform)
W3 <- data.frame(mag = 3, female = 1, magAnalysisWoman = 3 * 1)
## Scenarios: Men, M = 2
M2 <- data.frame(mag = 2, female = 0, magAnalysisWoman = 2 * 0)
## Scenarios: Men, M = 3 (post-reform)
M3 <- data.frame(mag = 3, female = 0, magAnalysisWoman = 3 * 0)
# Women and Men Data
Wdata <- rbind(W2, W3)
Mdata <- rbind(M2, M3)

## Simulate betas and vcov
set.seed(123)
sims1 <- rmvnorm(n = 1000, mean = coef(m1), sigma = vcov1)
sims2 <- rmvnorm(n = 1000, mean = coef(m2), sigma = vcov2)
sims3 <- rmvnorm(n = 1000, mean = coef(m3), sigma = vcov3)
sims4 <- rmvnorm(n = 1000, mean = coef(m4), sigma = vcov4)
sims5 <- rmvnorm(n = 1000, mean = coef(m5), sigma = vcov5)
sims6 <- rmvnorm(n = 1000, mean = coef(m6), sigma = vcov6)
sims7 <- rmvnorm(n = 1000, mean = coef(m7), sigma = vcov7)
sims8 <- rmvnorm(n = 1000, mean = coef(m8), sigma = vcov8)

## Calculate Predicted Values
# Women
pred1_W <- apply((sims1 %*% t(Wdata[2,])) - (sims1 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred2_W <- apply((sims2 %*% t(Wdata[2,])) - (sims2 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred3_W <- apply((sims3 %*% t(Wdata[2,])) - (sims3 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred4_W <- apply((sims4 %*% t(Wdata[2, -2])) - (sims4 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) # omit female
pred5_W <- apply((sims5 %*% t(Wdata[2, -2])) - (sims5 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) # omit female
pred6_W <- apply((sims6 %*% t(Wdata[2, -2])) - (sims6 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) # omit female
pred7_W <- apply((sims7 %*% t(Wdata[2, -2])) - (sims7 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) # omit female
pred8_W <- apply((sims8 %*% t(cbind(1, Wdata)[2,])) - (sims8 %*% t(cbind(1, Wdata)[1,])), 2, quantile, c(0.975, 0.5, 0.025)) # add intercept

# Men
pred1_M <- apply((sims1 %*% t(Mdata[2,])) - (sims1 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred2_M <- apply((sims2 %*% t(Mdata[2,])) - (sims2 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred3_M <- apply((sims3 %*% t(Mdata[2,])) - (sims3 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred4_M <- apply((sims4 %*% t(Mdata[2, -2])) - (sims4 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) # omit female
pred5_M <- apply((sims5 %*% t(Mdata[2, -2])) - (sims5 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) # omit female
pred6_M <- apply((sims6 %*% t(Mdata[2, -2])) - (sims6 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) # omit female
pred7_M <- apply((sims7 %*% t(Mdata[2, -2])) - (sims7 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) # omit female
pred8_M <- apply((sims8 %*% t(cbind(1, Mdata)[2,])) - (sims8 %*% t(cbind(1, Mdata)[1,])), 2, quantile, c(0.975, 0.5, 0.025)) # add intercept

# Building blocks of table C.10
bill_hat_women <- c(round(pred1_W[2,], 2), round(pred2_W[2,], 2), round(pred3_W[2,], 2),
	round(pred4_W[2,], 2), round(pred5_W[2,], 2), round(pred6_W[2,], 2),
	round(pred7_W[2,], 2), round(pred8_W[2,], 2))
ci_women <- c(paste0('(', round(pred1_W[3,], 2), ', ', round(pred1_W[1,], 2), ')'), 
	paste0('(', round(pred2_W[3,], 2), ', ', round(pred2_W[1,], 2), ')'), 
	paste0('(', round(pred3_W[3,], 2), ', ', round(pred3_W[1,], 2), ')'), 
	paste0('(', round(pred4_W[3,], 2), ', ', round(pred4_W[1,], 2), ')'), 
	paste0('(', round(pred5_W[3,], 2), ', ', round(pred5_W[1,], 2), ')'), 
	paste0('(', round(pred6_W[3,], 2), ', ', round(pred6_W[1,], 2), ')'), 
	paste0('(', round(pred7_W[3,], 2), ', ', round(pred7_W[1,], 2), ')'), 
	paste0('(', round(pred8_W[3,], 2), ', ', round(pred8_W[1,], 2), ')'))


bill_hat_men <- c(round(pred1_M[2,], 2), round(pred2_M[2,], 2), round(pred3_M[2,], 2),
	round(pred4_M[2,], 2), round(pred5_M[2,], 2), round(pred6_M[2,], 2),
	round(pred7_M[2,], 2), round(pred8_M[2,], 2))
ci_men <- c(paste0('(', round(pred1_M[3,], 2), ', ', round(pred1_M[1,], 2), ')'), 
	paste0('(', round(pred2_M[3,], 2), ', ', round(pred2_M[1,], 2), ')'), 
	paste0('(', round(pred3_M[3,], 2), ', ', round(pred3_M[1,], 2), ')'), 
	paste0('(', round(pred4_M[3,], 2), ', ', round(pred4_M[1,], 2), ')'), 
	paste0('(', round(pred5_M[3,], 2), ', ', round(pred5_M[1,], 2), ')'), 
	paste0('(', round(pred6_M[3,], 2), ', ', round(pred6_M[1,], 2), ')'), 
	paste0('(', round(pred7_M[3,], 2), ', ', round(pred7_M[1,], 2), ')'), 
	paste0('(', round(pred8_M[3,], 2), ', ', round(pred8_M[1,], 2), ')'))

tableC10 <- rbind("Women "= bill_hat_women, " " = ci_women, 
	"Men "= bill_hat_men, " " = ci_men)
colnames(table2) <- paste0('Model ', 1:8)
stargazer(table2,
	title = "Predicted percentage point change in the cosponsorship portfolio on women’s issues when M increases by one")
